A Gradient Microstructure Improves the Barrier Properties of Flake-Filled Composite Films: A Computational Study

Composite films of a graded miscrostructure hold the promise of achieving optimal use of the filler material, resulting in composites with improved and application-taylored properties. In the context of barrier materials in which the reinforcing phase comes in the form of flakes or platellets, concentrating the filler particles in certain critical regions is thought to achieve economy in filler usage while ensuring superior barrier performance. The objective of the present article is to quantitatively test this hypothesis and provide guidelines on the expected barrier improvement. A model is developed, according to which a graded miscostructure in a composite film offers a quantitative improvement over an equivalent homogeneous microstructure; this improvement is quantified using a coefficient β, which depends on the form of the graded miscrostructure, specifically the distribution of the number-density of the filler particles across the film. It is shown that β=1 for a uniform microstructure and β>1 for a graded one, indicating that a graded miscrostructure will indeed result in improved barrier properties. Analytical expressions for β are developed for certain typical distributions; for a linear filler distribution, it is shown that β=4/3. This model is tested against detailed multi-particle simulations and is found to be in excellent agreement with computational results.


Introduction
Materials of a gradient microstructure, otherwise known as functionally gradient materials (FGMs), form a class of heterogeneous materials that exhibit a gradual transition of composition or morphology between two extreme states [1]. One example of compositional grading is seen in metal-matrix composite FMGs used in coatings or components designed to afford thermal or wear protection as well as minimize the effect of thermal stresses or the propagation of cracks; in these, a ceramic-rich exterior is seamlessly connected to a metallic interior by a material having a graded composition [2]. An example of microstructural grading can be found in the case of articular cartilage, a specialized connective tissue found in synovial joints. In this, the orientation of the constituent collagen fibers changes gradually from parallel to perpendicular to the articular surface [3,4], resulting in an optimal distribution of stress. FGMs abound in nature, as, for example, in the structures forming human bone, in several sea-shell tissues as well as in plants. The earliest reference to manufactured FGMs appeared in the late 1980s with the pioneering work of Japanese researchers [5,6] in the context of the development of thermal protection materials for aerospace applications. Since then, many applications of FGMs have appeared in bio-medical devices (prosthetics, for example, in the form of various HAP/Ti composite coatings [7][8][9]), nuclear energy, and other high-temperature applications [10,11], sensors and actuators, electricity generation, optical components, and many others. The advent of additive manufacturing and the flexibility it offers to gradually vary the composition or morphology in a component at the microstructure level has resulted in a renewed interest in FMGs in recent years [12], while a recent review of FGMs as related to manufacturing and properties can be found in [13]. In the area of micro-scale heat transfer, the effect of a graded microstructure has been studied, computationally and experimentally, by Dondero and co-workers [14]; the methodology, both computationally (the use of the fast multiple boundary element method, e.g., [15]) and experimentally (holes drilled in an aluminum plate, thus obtaining a variety of gradient profiles) can be extended to construct and analyze more general model systems of the gradient microstructure.
Contrary to the above applications, the concept of a gradient structure has not yet been fully utilized in the area of flake composites, even though multi-layer applications involving flake composites appear, for example, in sound insulation using mica flakes [16]. Flake composites, in which the reinforcing particles have the form of platelets (flakes), are characterized by a high aspect ratio and the very high area-to-volume ratio of the reinforcing particles. In mass-transfer applications, these particles act as local impermeable barriers and, therefore, force diffusing solute molecules to follow a highly tortuous path as they diffuse through the material [17]. This results in significantly reduced effective diffusion coefficients, and thus flake composites have found usage in numerous applications in which the diffusion rate of gasses through a film needs to be minimized, such as in corrosion [18], fire resistance [19], and packaging [20]. The prediction of the effective diffusion coefficient D e f f in such composites is well advanced [21][22][23][24][25][26][27][28][29], as in the case when the flakes are uniformly distributed in the composite . However, the potential effect of a graded distribution of the flakes on their barrier effectiveness has not been studied. The objective of this study is to fill this gap and provide estimates for the barrier improvement factor (BIF = D 0 /D e f f ) of gradient flake composites for the base-line case of flakes aligned normally to the direction of macroscopic diffusion. To our knowledge, this is the first study that deals with flake composites of a graded structure. The paper is organized as follows: In Section 2, we develop a novel closed-form model for the BIF of composites having various common graded structures. In Section 3, a detailed computational analysis of the relevant systems is presented. Section 4 summarizes the computational results and compares them to the predictions of the analytical model.

A Model for the Effective Barrier Improvement Factor (BIF) of a Graded Composite Film
Consider a distribution of flakes, as shown in the schematic of Figure 1, in which N(y) is the local density of the flakes at a location y of height δy. This is related to the total number of flakes across the unit cell through N T = H 0 N(y)dy. At a location y, the local concentration measure will be αφ(y), where φ is the local volume fraction of flakes in a horizontal cross section and α is the flake aspect ratio, defined as the ratio between flake length l and flake thickness t. It can be shown [29] that the product αφ can be expressed as αφ(y) = (l/H) 2 N(y). Using this, the local BIF, if expressed by the model of Lappe et al. [25], will be where D 0 is the diffusion coefficient of the solute molecules in the matrix material; D e f f is the effective diffusion coefficient of the same molecules in the composite; and λ is a geometrical factor, expressing the (increased) path a solute molecule is forced to take as it navigates its way around a flake. The model of Lappe at al. [25] has been found to be in excellent agreement with the results of 2D and 3D computations in spatially uniform composites and, for 2D systems, the geometrical constant λ has been found to take the value of λ = 2.5 [30,31]. We will use this value in this study. In light of Equation (1), an average value of the BIF across the film will be The geometry is periodic in the horizontal direction, and periodic boundary conditions are applied on the vertical boundaries of the unit cell.
In this example, αφ = 0.5 where αφ is the macroscopic average flake concentration. In total, 500 flakes were placed in the unit cell so that their vertical positions followed a normal distribution with maximum density at 0.5H and σ = 0.2.
Equation (2) is based on the assumption that mass transfer across the gradient composite can be described by a resistance-in-series model. This is a correct assumption for layered composites containing randomly placed flakes, in which no low resistance paths form through flake layers. Expanding the right-hand-side of Equation (2), we obtain Obviously, N T = αφ(H/l) 2 and therefore the second term in the right-hand side of Equation (3) is (2/λ)(αφ). The third term in the right hand side of Equation (3) can be written as where is a positive factor that depends only on the form of the distribution N(y); evidently, β expresses an improvement in barrier effectiveness brought about by the graded distribution of the flakes. Through Equations (3) and (4), the effective BIF can be expressed as According to Equation (6), the BIF of a graded composite film is expected to be higher than the BIF of an equivalent (same number and dimensions of flakes per unit area/volume) composite in which there is no graded structure, but, instead, the flakes are uniformly distributed throughout the film thickness. Obviously, in the latter case, it can easily be shown that αφ(y) = αφ and β = 1. Notice that, from Equation (6), it is clear that β operates on the quadratic term of the BIF model; therefore, its effect will be seen best in concentrated systems. Conversely, grading is only expected to have a minor effect in dilute flake dispersions since, in such systems, BIF is dominated by the linear term in (αφ).
Equation (6) is significant since, if verified, it provides a direct means of assessing the increase in barrier effectiveness caused by grading and helps to identify microstructures that offer the best barrier enhancement. In the following, we will perform computer simulations in order to determine the range of validity of Equation (6). We will consider the following four types of flake distributions: (i) A half-normal distribution, according to which the flake concentration is maximum on the surface (y = 1) and minimum at the bottom (y = 0) of the film.
In this case, the number-density of the flakes N(y) is given by where the constant A is related to the total number of flakes N T through a total balance: (5), it can be shown that in this case the enhancement factor β is a function of the steepness of the distribution (expressed by the parameter σ). The relation between β and σ is: It can be formally shown that as σ → ∞, then β → 1. This result has also been verified with numerical computations, as will be shown later in this paper. (ii) A linear distribution between the two sides of the film. In this case, it is N(y) = ky, where k is a constant. Since N T = H 0 kydy, it follows that k = 2N T /H 2 , and therefore N(y) = (2N T /H 2 ) × y. Taking into account that αφ(y) = (l/H) 2 N(y) and introducing this in Equation (5), we derive (iii) A triangular distribution in which the flake concentration is maximum at the center and drops linearly towards the two edges of the film. It can be shown that in this case it is also β = 4/3. (iv) A normal distribution, with the maximum flake density occurring at the center of the film (y = 0.5H). In this case, the number-density of the flakes varies as In a manner similar to the one used in the case of the half-normal distribution, the enhancement factor β can be shown to depend only on the variance (σ) of the distribution; the result is: It can be shown that as the variance of the distribution increases, the enhancement factor β approaches 1 asymptotically .

Computational
Computations are carried out in two-dimensional unit cells of square shape (side length H), which are formed by placing N T flakes according to a desired distribution N(y). The flakes are unidirectional, with their plane perpendicular to the direction of the macroscopic diffusion (y). We are interested in the case of thin flakes, whose crosssections appear as 1D lines, parallel to each other and perpendicular to the direction of the macroscopic diffusion, as shown in Figure 2. When the flakes are aligned with their axes perpendicular to the direction of diffusion, the hindrance to diffusion (and thus the BIF) are maximized. For this reason, this form of flake alignment is the preferred alignment in flake composites and is the configuration mostly studied, computationally as well as theoretically. Most resin processing/application methods for polymer melts, epoxy resins, or paints involve flows with substantial elongational or shear components. These flows are known to align the flakes in the direction of flow [32], which in practical applications is the same as the direction of the coated substrate and normal to the direction of macroscopic diffusion.
The state of a flake dispersion is characterized as dilute, semi-concentrated, or concentrated, based on the value of the product αφ, where α = l/t is the flake aspect ratio and φ is the flake volume fraction. If αφ > 1, we are in the concentrated regime. If 1/α < αφ < 1, we are in the semi-concentrated regime, while a dispersion in which αφ < 1/α is considered dilute. Additionally, since in our case the placement of flakes is non-uniform, the local αφ varies, and at various positions it can exceed the above characterization limits. This is especially evident at smaller values of σ since these lead to more dense areas inside the unit cell due to higher concentration of flakes. So, even if a film is characterized by an average value of αφ, locally it shows a different character, and this affects the overall barrier performance. This is ilustrated in Figure 3.
As shown in Figure 4, the procedure that is used for running simulations for each set of geometrical parameters consists of the following steps:

1.
A set of parameters is determined, and the in-house application generates the appropriate geometry files. A unit cell is created with the desired dimensions, and for each set of φ and distribution the number of flakes is calculated. Then, the flakes are placed using a random sequential addition algorithm [33], according to which random numbers define the coordinates of the center of each line. The RSA algorithm performs checks to ensure no flake overlap occurs, and to ensure that a minimum separation distance between neighboring flakes is maintained. The latter is essential so that a computational mesh can be generated in the narrow gaps separating adjacent flakes. After the completion of the unit cell calculation, the case directories are created with the appropriate scripts for execution. A master script is created that will run in the main node and coordinate the execution at the various nodes.

2.
The master script is executed; it uses a greedy strategy to assign the cases to each computational node. According to this strategy, each node in the computational cluster sequentially scans the case directories until an 'unsolved' case is found. In this case, the actions in step 3 are executed. If all of the cases are scanned and no other 'unsolved' cases are found, the script terminates its execution.

3.
The nodes are searching for the existence of status files inside each directory that indicate the state of the case. If the status is 'unsolved', it is flagged as 'running' and the node starts the execution of a sub-script that runs the simulation. If a different node accesses that case folder, it will read the status and will skip it if its status is any other than 'unsolved'. If a case is 'unsolved', then the actions of step 4 are executed.

4.
A sub script is run in each case directory that executes the commands for running the simulation. The first step is the creation of the mesh. Then the mesh is converted using OpenFOAM utilities, and the boundary conditions are applied. The flakes are being simulated as being one dimensional entities that have only length or 'baffles', as referenced in the boundary conditions settings. In the case of baffles, boundary conditions are applied to the internal faces of both cells that lie on each side of the geometrical entity, effectively creating an internal boundary inside the mesh. The simulation is split in order to run in parallel using the the public domain openMPI implementation of the standard message passing interface (MPI). After a solution is achieved, the results are written in a file and the case is flagged as 'solved'. 5.
The in-house application is used to scan the case directories and is retrieving the simulation results. Here, it is interesting to note that there is no need to finish all of the simulations before the results are obtained. The results that have been obtained at any time can be collected as needed.  In step #2, the case files are uploaded to the cluster and the master script is executed at a master node sending the individual scripts to the computational nodes. Each node runs the cases in step #3 using a greedy algorithm. Every script that runs (step #4) creates the geometry using boundary representation geometry and then creates the mesh. Using OpenFOAM utilities, the mesh is converted, boundary conditions are applied, and the discretized model is partitioned according to the physical cpu cores in every node and executed in parallel. Finally, the results are obtained using again the OpenFOAM utilities. At step #5, the results are collected using in-house utilities for further analysis.
These simulations were run in the Linux environment of the National HPC facility ARIS using 8 nodes with 16 cpus each.
The unit cells are periodic in the x-direction. The boundary conditions for concentration are cyclic in the x-direction, namely, C(0, y) = C(L, y). In the y-direction, a fixed concentration gradient ∆C is imposed; therefore, fixed concentration values are prescribed on the top and bottom boundaries, simulating the steady-state case of a film placed as a barrier between two regions. The flakes are considered to be impermeable, and thus at every point on their surface we have dC/dn = 0. At each level of αφ, we generate 10 different geometries, which differ from one another only by the spatial placement of the flakes. The files for each case setup and the geometric information regarding the unit cells is created using in-house developed applications after setting the flake distribution characteristics and the rest of the case settings. The geometric files are imported in the meshgeneration application GMSH [34]; the generated meshes consist of around ∼10 7 triangular elements. More information on mesh generation and spatial mesh convergence is included in the Appendix A. The case files created above are then imported into OpenFOAM [35] and the Laplace equation, namely, is solved using the well-tested solver LaplacianFoam. The adoption of Equation (12) implies that the matrix material is assumed to be homogeneous throughout, also including the regions adjacent to the flakes. This is a simplification for certain material/flake combinations as it is known that the presence of flakes may alter the conformation/crystallization and thus the properties of the matrix in their vicinity [36,37]. The solution to Equation (12) gives the value of dC/dn at each point on the top/bottom boundaries, from which the total mass flux is computed as Knowledge of the mass flux allows us to compute an effective diffusion coefficient of the composite film as Example concentration profiles are shown in Figure 5.

Half-Normal Distribution of Flakes
Equation (6) suggests that if the computed values of the BIF of the (graded) composite film are plotted as (BIF − 1)/αφ vs. αφ, they should fall on a straight line, whose slope will be equal to β/λ 2 . This is indeed the case, as demonstrated in Figure 6, which shows that the enhancement factor (β) increases with increased grading (decreasing value of σ), and this increase is more evident for larger values of αφ. The numerical results for β, obtained from Figure 6, are plotted in Figure 7 and also compared to the predictions of Equation (8).  It is clear that β asymptotically approaches the value of 1 with increasing σ − as expected since for σ → ∞, the composite approaches a homogeneous state. It is also clear that the numerical results for β and the predictions of Equation (8) deviate at lower values of σ, that is, for films showing very high grading. The reason for this discrepancy can be found in the assumption underlying the derivation of Equation (6), namely, that the average BIF can be computed as a sum of resistances in a series. This assumption implies that there are no preferential (low resistance) paths in the flake layers. Such paths, however, can sometimes be found in very steep flake concentration profiles, in which case the great majority of the flakes are concentrated in a very thin subsurface layer, as can be seen in Figure 8.

Linear Distribution
Equation (9) predicts that in the case of a linear variation of flake number-density, the enhancement factor is β = 4/3. In this case, plotting the results for the BIF as in the previous section, namely, (BIF −1 )/αφ vs. αφ, should result in one straight line, whose slope should be β/λ 2 = 0.21-if we take λ = 2.5. Figure 9 shows the results, which form a straight line with slope 0.209, resulting in a value of β equal to 1.306, compared to a value of 4/3=1.333 expected from theory.

Triangular Linear Distribution
As shown earlier, the value of β in this case should be 4/3. The computational results are shown in Figure 10. In this case, the slope of the regression line is 0.207. It is interesting to note here the similarities between the linear and triangular distribution. From the computational results, we see that they are virtually identical. This was to be expected, since these two distributions are geometrically similar. The linear distribution can be considered as a overlapping of two areas with triangular distributions. This is demonstrated in Figure 11, where we can see that a unit cell with linear distribution can be produced from a unit cell with triangular distribution.  Figure 12 is a summary of the results for the BIF in the case of a normal distribution of flakes, with the maximum concentration located at the center of the film as a function of αφ and σ. It can be seen that increasing σ beyond 0.3 results in very similar BIF vs. αφ curves; this is an indication that the gradient structure of the film will have an appreciable effect on its barrier properties only for relatively steep concentration gradients, typically σ < 0.2. Figure 12 also provides a guide for film design; if a flake-filled film with a target BIF is to be obtained, this can be achieved using various total flake amounts, depending on the extent of grading. In the example of Figure 12, a BIF of 20 can be achieved with αφ = 4.3 if σ = 0.05 (an extreme case of grading); this becomes 5.8 at σ = 0.1 and 6.8 at σ = 0.15. To achieve a BIF = 20 with uniformly distributed flakes, the required value of αφ would have to be λ( √ 20 − 1) = 8.7, which is double the flake amount if grading with σ = 0.05 was used and 50% higher if a grading with σ = 0.1 was used. The enhancement factor β was predicted to be a function of σ, as expressed by Equation (11). Figure 13 shows a comparison between the predictions of Equation (11) and the computational results for β.

Conclusions
Theoretical analysis and detailed computational results point to the conclusion that the barrier factor of a flake-filled film will increase if the spatial distribution of the flakes is non-uniform in the direction of macroscopic diffusion. The extent of improvement is directly linked to the form of the distribution through a coefficient (β) we introduce in this study. This coefficient operates on the quadratic term of the BIF model; therefore, the impact of grading is expected to be more significant in the semi-dilute and concentrated regimes. Analytical expressions for this coefficient are obtained for the cases of typical distributions such as linear, half-normal, and normal. From computational results, the coefficient β is also found and the obtained values are in very good agreement with the analytical expressions. Computational results also confirm that the highest impact of grading will be observed in concentrated films. Finally, our results are obtained for the case of fully-aligned flakes; for composites in which the flakes assume random orientations or in which the flakes are only partially aligned, the BIF reported here will form an upper bound [38].
The values of BIF for different flake distributions and parameters are shown in Figure 14. We can see that the BIF values obtained from a half normal distribution with small σ are the highest. This is explained considering that the half normal distribution becomes very dense in the upper surface of the unit cell. This is also true for the normal distribution with σ = 0.05 since this distribution also has a very dense area but this time in the middle of the unit cell. The triangular distribution and half triangular distribution (not shown in Figure 14 but giving almost identical results with the triangular distribution as shown in Figure 10) give the smallest increase in the BIF.
From the above it is concluded that in order to achieve high values of BIF, a manufacturing process which creates materials with highly flake-dense zones should be followed. If this is not achieved, then the BIF will effectively be the same as that obtained from a material with a random distribution of flakes. Small deviations from maximum grading can lead to a large reduction in BIF, as can be seen in the normal distribution where a difference from σ = 0.05 to σ = 0.15 can result in a two-fold reduction in BIF, making the barrier properties of the material similar to those of a material having a random distribution of flakes.

Conflicts of Interest:
The authors declare no conflict of interest.
An example of the application of the REM in our problem is given in the following, for a case with αφ = 1.0. For a given geometry, three different meshes are made using a refinement factor r = 1.8. The results of the simulations are summarized in Figure A1.
In this typical case, the order of convergence is found using Equation (A3) and is equal to 3.6226. Using Equation (A4), the effective diffusivity (D e f f ) is 0.49826. Since we used three grids with an effective ratio r = 1.8, the S f is equal to 1.25. The ratio was kept constant by utilizing the setting of 'CharacteristicMeshLength' of GMSH. In order to calculate the asymptotic rate of convergence, we use Equations (A5) and (A6), and the resulting values are B 1,2 = 0.08168 and B 2,3 = 0.68358. The ratio B 2,3 /B 1,2 = 0.99518 ∼1. Therefore, from the above analysis and comparing the results from the simulation with >10 7 grid elements, we can see that we are well within the asymptotic range. By performing a similar analysis for the number of flakes, we end up using N = 4000 for our simulations. . REM was carried out for all simulation settings, and as a result it was concluded that meshes with around 10 7 elements and flake number N = 4000 are sufficient for all αφ and σ levels examined. (b) Meshes created from the above procedure from coarse mesh (normalized grid spacing ∼10) to finer mesh (normalized grid spacing ∼1). A local size refinement was applied at the flake edges to more efficiently resolve the concentration gradients in these areas.